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VARIATIONAL MULTISCALE PROPER ORTHOGONAL 
DECOMPOSITION: CONVECTION-DOMINATED 
CONVECTION-DIFFUSION EQUATIONS 

TRAIAN ILIESCU AND ZHU WANG 

Abstract. We introduce a variational multiscale closure modeling strategy 
for the numerical stabilization of proper orthogonal decomposition reduced- 
order models of convection-dominated equations. As a first step, the new 
model is analyzed and tested for convection-dominated convection-diffusion 
equations. The numerical analysis of the finite element discretization of the 
model is presented. Numerical tests show the increased numerical accuracy 
over the standard reduced-order model and illustrate the theoretical conver- 
gence rates. 



1. Introduction 

One of the most successful dynamical systems ideas in the study of turbulent 
flows has been the Proper Orthogonal Decomposition (POD) |16[ I37j . POD has 
been used to generate reduced- order models (ROMs) for the prediction and control 
of structure dominated turbulent flows [TJ [U [5J [321 133] • The idea is straightforward: 
Instead of using billions of local finite element basis functions equally distributed 
in space, POD uses only a few (usually O(10)) global basis functions that represent 
the most energetic structures in the system. Thus, the computational cost in a 
direct numerical simulation (DNS) of a complex flow can be reduced by orders of 
magnitude when POD is employed. 

Despite its widespread use (hundreds of papers being published every year), POD 
has several well-documented drawbacks. In this report, we address one of them, 
namely the numerical instability of a straightforward POD Galerkin procedure ap- 
plied to a complex flow 2J. To address this issue, we draw inspiration from the 
methodologies used in the numerical stabilization of finite element discretization 
of convection-dominated flows. Specifically, we employ the variational multiscale 
(VMS) approach used by Layton in [25], which adds artificial viscosity only to the 
smallest resolved scales. We also note that an approach similar to that used in [25] 
was proposed by Guermond in [13l [14] . 
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We emphasize that the VMS philosophy is particularly appropriate to the POD 
setting, in which the hierarchy of small and large structures appears naturally. 
Indeed, the POD modes are listed in decreasing order of their kinetic energy content. 
We also note that, although a VMS-POD approach was announced in [51 [7] and 
another VMS-POD approach was used in [4], to the authors' knowledge this is the 
first time that the VMS formulation in 29 has been applied in a POD setting. 

In this report, the new VMS-POD model is analyzed and tested in the numerical 
approximation of a convection- dominated convection- diffusion problem 

!u t - e Au + b • Vu + gu = f in (0, T] X ft , 

u(x,0) = u (x) in ft, 

u(x,t) = on (0, T] x <9ft , 

where e <C 1 is the diffusion parameter, b with ||b|| = 0(1) the given convective 
field, g the reaction coefficient, / the forcing term, ft C R 2 the computational do- 
main, t G [0, T], with T the final time, and ito(") the initial condition. Without loss 
of generality, we assume in what follows that the boundary conditions are homo- 
geneous Dirichlet. We emphasize that the new VMS-POD model targets turbulent 
flows described by the Navier- Stokes equations (NSE). We chose the mathematical 
setting in ( |1.1[ ), however, because it is simple, yet relevant to our ultimate goal 
(since e <C ||b||). Of course, once we fully understand the behavior of the new 
VMS-POD model in this simplified setting, we will analyze and apply it in the 
NSE setting. 

The rest of the paper is organized as follows. In Section[2j we briefly describe the 
POD methodology and introduce the new VMS-POD model. The error analysis 
for the finite element discretization of the new model is presented in Section [3] The 
new methodology is tested numerically in Section [4] for a problem displaying shock- 
like phenomena. Finally, Section [5] presents the conclusions and future research 
directions. 



2. Variational Multiscale Proper Orthogonal Decomposition 

2.1. Proper Orthogonal Decomposition. In this section, we briefly describe 
the POD. For a detailed presentation, the reader is referred to [TBI |2"TI 157] . 

Let X be a real Hilbert space endowed with the inner product (•, -)x, and tt(-, t) G 
X, t G [0, T] be the state variable of a dynamical system. Given the time instances 
t\, . . . , tjv G [0, T], we consider the ensemble of snapshots 



(2.1) 



V := span{u(-,ti), . . .,u(-,t N )} . 



with dim V = d. The POD seeks a low-dimensional basis {<fi, ■ ■ ■ , 'fir}, with r <C rf, 
which optimally approximates the input collection. Specifically, the POD basis 
satisfies 



(2.2) 



1 N 

-Y 



X 



subject to the conditions that (tfi,ipj)x = 1 < i,j < r. In order to solve (2.2), 
we consider the eigenvalue problem 



(2.3) 



K v = A v , 
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where K € R 



NxN 



1 



with Kij — — (u(-,tj),u(-,ti)) x , is the snapshot correlation 
matrix, Ai > A2 > . . . > A<j > are the positive eigenvalues, and v^, k = 1, . . . , d 
are the associated eigenvectors. It can then be shown (see, e.g., [SI HZ])) that the 
solution of (2.2 1 is given by 

1 N 

(2-4) ¥>*(•) = -7i= J>fc)i ±<k<r, 



where (nfe)j is the j-th component of the eigenvector Vk- It can also be shown that 
the following error formula holds: 

2 

I N r 

i=l 



(2.5) 



N 

«=i 



A" 



d 

E 

j=r+l 



In what follows, we will use the notation X r — span{(^!, ip 2 , ■ ■ ■ , <p r } ■ Although X 
can be any real Hilbert space, in what follows we consider X := Hq(£1). 

In the form it has been presented so far, POD seems to be only a data com- 



pression technique. Indeed, equation (2.2 1 simply says that the POD basis is the 
best possible approximation of order r of the given data set. In order to make 
POD a predictive tool, one couples the POD with the Galerkin procedure. This, in 
turn, yields a ROM, i.e., a dynamical system that represents the evolution in time 
of the Galerkin truncation. We now briefly present the derivation of this ROM, 
we highlight one of its main drawbacks and we propose a method to address this 
deficiency. 

The POD-Galerkin truncation is the approximation u r <E X r of u: 



(2.6) 



u r (x, t) 



3=1 



;(*) Pi 0» 



Plugging (2.6) into (1.1) and multiplying by test functions in X r C X yields the 



POD-Galerkin (POD-G) model 

(2.7) (ur,t, v r) + e(Vw r , Vv r ) + (b • V« r , v r ) + (g u r , v r ) — (/, v r ) Vw r el' 



The main advantage of the POD-G model (2.7) over a straightforward finite element 



(FE) discretization of ( 1.1 ) is clear - the computational cost of the former is dramat- 



ically lower than that of the latter. There are, however, several well-documented 



disadvantages of (2.7), such as its numerical instability in convection-dominated 
flows |36) . To address this issue, we draw inspiration from the methodologies used 
in numerical stabilization of finite element discretizations of such flows. 

2.2. Variational Multiscale. The Variational Multiscale (VMS) method intro- 
duced by Hughes and his group [HI HH [20l [21] has been successful in the nu- 
merical stabilization of turbulent flows pH H2 [22l [23l EH ES] . The idea in VMS 
is straightforward: Instead of adding artificial viscosity to all resolved scales, in 
VMS artificial viscosity is only added to the smallest resolved scales. Thus, the 
small scale oscillations are eliminated without polluting the large scale components 
of the approximation. The VMS method has been extensively developed, various 
numerical methods being used. The finite element discretization of the resulting 
VMS model has evolved in several directions: Hughes and his group proposed a 
VMS formulation for the NSE in which a Smagorinsky model [5[ [35] was added 
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only to the smallest resolved scales [HI [HI [2Ql EQ ■ A different type of VMS ap- 
proach, based on the residual of the NSE, was proposed in Bazilevs et al. [31. One 
of the earliest VMS ideas for convection-dominated convection-diffusion equations 
was proposed by Guermond in [HI [T3]. In this VMS formulation, the smallest 
scales were modeled by using finite element spaces enriched with bubble functions. 
Layton proposed in [29] a VMS approach similar to that of Guermond. In this VMS 
approach, however, the smallest resolved scales were modeled by projection on a 
coarser mesh. The VMS approach proposed in [5S] was extended to the NSE in a 
sequence of papers by John and Kaya [HI [Ml 122] • The variational formulation 
used by the FE methodology fits very well with the VMS approach. The definition 
of the smallest resolved scales, however, often poses many challenges to the FE 
method. Indeed, one needs to enrich the FE spaces with bubble functions [T31 114] . 
consider hierarchical FE bases [TH], or use a projection on a coarser mesh [25] , 

2.3. The VMS-POD Model. POD represents the perfect setting for the VMS 
methodology, since the hierarchy of the basis is already present. Indeed, the POD 
basis functions are already listed in descending order of their kinetic energy con- 
tent. Based on this observation, we next propose a VMS based POD model. To 
this end, we consider the following spaces: X :— Hq(CI), X h c Hq(CI), X r := 
span{(^i, tp 2l ■ ■ ■ , (fir}, X R :— span{</?i, tpz, . . . , <pn}, where R < r, and L , where 
L R will be defined later. Note that X R C X r C X h c X. We also consider 
Pr : L 2 (0) — > L R , the orthogonal projection of i 2 (il) on L R , defined by 

(2.8) (u-P R u,v R )=Q Vv R eL R . 

Let also P R := I — P R . We are now ready to define the Variational Multiscale 
Proper Orthogonal Decomposition (VMS-POD) model: 

^ ^ (u r ,t, IV ) + £ (V«r, Vw r ) + a (P R Vu r , P R Vv r ) 

+ (b • Vu r , v r ) + (gu r , v r ) — (/, v r ) Vtv € X r . 



The third term on the LHS of ( |2.9[ ) represents the artificial viscosity that is added 
only to the smallest resolved scales of the gradient. We note that, although a 
VMS-POD approach was announced in [51 [7] and another one was used in [3] , to 
the authors' knowledge this is the first time that the VMS formulation in 29J is 
applied in a POD setting. 

In the next two sections, we will first estimate the error made in the finite clement 



discretization of the new VMS-POD model (2.9 1 and then use it in a numerical test 



3. Error Estimates 

1 N 

In this section, we prove estimates for the average error — — - | 

' n=0 



where the approximation u n is the solution of (3.11) (the weak form of ( |1.1[ )) and 



< is the solution of ( |3.12[ ) (the FE discretization of the VMS-POD model p^9|). 
To this end, we follow the approach in [TS] (see also [35]). We emphasize, however, 
that our presentation is different in that it has to include several results pertaining 
to the POD setting. To this end, we use some of the developments in [3T] (see also 

[UlEHlEolEl]). 

We start by introducing some notation and we list several results that will be 
used throughout this section. For clarity of notation, we will denote by C a generic 
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constant that can depend on all the parameters in the system, except on d (the 
number of POD modes retained in the Galerkin truncation), N (the number of 



snapshots), r (the number of POD modes used in the POD-G model (2.7)), R (the 



number of POD modes used in the projection operator in the VMS-POD model 



(2.9)), h (the mesh-size in the FE discretization), a (the artificial viscosity coeffi- 
cient), and e (the diffusion coefficient). Of particular interest is the independence 
of the generic constant C from e. Indeed, we will prove estimates that are uniform 
with respect to e, which is important when convection-dominated flows (such as the 
NSE) are considered. 

We introduce the bilinear forms b(u, v) := (b • Vw,u) + (gu,v), a(u,v) := 
e(Vu, Vv)+b(u, v), and A(u, v) :— a(u,v) + a (P R \7u, P R \7v). We also consider the 
weighted norm b a :— a \\u\\ 2 + b || Vu|| 2 +a \\P R Wu\\ 2 . We now make the follow- 
ing assumption, which is used in proving the well-posedness of the weak formulation 
of ( fl~T| ). 

Assumption 1 (Coercivity and Continuity). 

(3.1) 5-^V-b>/3>0 and max{||g||, ||b||} = 7 > 0. 



For the FE discretization of (1.1), we consider a family of finite dimensional 
subspaces X h of X = H^(Q,), such that, for all v € H m+1 n X, the following 
assumption is satisfied. 

Assumption 2 ( Approximability) . 

(3.2) inf {\\v-v h \\+h\\Wv-Wv h \\} < C h m+1 \\v\\ m+1 1 < m < k, 
v h ex h 

where k is the order of accuracy of {X h }. We also assume that the finite element 
spaces {X h } satisfy the following inverse estimate. 

Assumption 3 (FE Inverse Estimate). 

(3.3) \\Vv h \\ < Ch- 1 \\v h \\ Vv h eX h . 



A similar inverse estimate for POD is proven in [28]. For completeness, we 
present it below. We also include a new estimate and present its proof. 

Lemma 3.1 (POD Inverse Estimate). Let M r S R rxr with Mij — (ipj,(fi) be the 
POD mass matrix, H r € R rxr with Hij = (V^j, V<£,) be the POD stiffness matrix, 
S r € M rxr with Sij — (ipj,(fii) H i be the POD mass matrix in the H x -norm, and 
|| • || 2 denote the matrix 2-norm. Then, for all v r G X r , the following estimates 
hold. 

(3-4) |k||z- < J\\Mr\\2\\Sr%\\Vr\\Hi, 

(3.5) ||» r || H i < y/^SM^%\\v r \\^, 

(3.6) IIVWrH^ < yJ\\H r \\2\\Mr% \\Vrh>- 
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Proof. The proof of estimates (3.4) and (3.5) was given in [28] (see Lemma 2 and 
Remark 2). The proof of (3.6 1 follows along the same lines: Let v r = YZ=i x jfi an d 
x = (afi, . . . ,x r ) T . From the definition of H r , it follows that ||Vu r ||| 2 = x T i? r x. 
Since H r is symmetric, its matrix 2- norm is equal to its Rayleigh quotient |10j : 
xT Thus, we get: 



H r 1 1 2 — max ■ 



(3.7) 



2 

r\\L 2 



ff r X < ||F r || 2 X T X. 



Furthermore, since M~ Y is also symmetric, we get y T M~ 1 y < ||M j r 1 || 2 y T y for 
all vectors y G R r . We also note that, since M r is symmetric positive definite, we 
can use its Cholesky decomposition M r = L r , where L r is a lower triangular 
nonsingular matrix [10] . Thus, letting y = L r x, we get: 



(3.8) 



\M-% 



> 



y T M-j 



X"^ L/^ {L r ) £ r ^ Z/ r x 
x T L T Lx 



T 
X X 

x T M r x 



Inequalities (3.7) and (3.8) imply the following inequality, which proves (3.6) 

IM" 1 



IVWrllt, < ||flr||2||M r - 



| 2 x T M r x = ||ff r || 2 



|2 



□ 



Remark 3.2. We note that, in our setting, (3.6) can be improved. Indeed, since S r 
is the identity matrix when X = Hq, we get: 



(3.9) \\Vv r \\ L 2 < \\v r \\ H i < VHSrllallAfr-^laHWrllia = VHM^IbllUrll^ 



We note, however, that in general (3.9) might not hold. 



To prove optimal error estimates in time, we follow 

u(t n )— u(t n _ i ) 



difference quotients du(t n ) = At , where n — 1, . . . , N, in the set of snap- 
shots V :— sp an \u (t Q ) 7 . . . , u(t N ),du(ti) 
error formula (2.5) becomes: 



and include the finite 
N, in the set 
, du(t?f)}. As pointed out in (35], the 



N 



(3.10) 



+ 



2N- 



1 



2N- 



j=0 



N 



1 ^ 

t=l 



u(-,ti) 



^2(u(;ti) , <Pj(-)) x <Pj(-) 
3=1 



X 



du(-,ti)-y2(du(-,ti), <Pj(')) x (pj(') 



3=1 



X 



= E ^ 

j=r+l 



After these preliminaries, we are ready to derive the error estimates. 

VveX. 



The weak form of (1.1) reads 
(3.11) 



(u t ,v) +a(u,v) = (f,v) 

The VMS-POD model for ( |3.11 1 with a backward Euler time discretization reads: 
Find u™ € V such that: 

(3.12) i«+ 1 -<,^) + A«+ 1 , Ur ) = (/»+ 1 )? ; r ) V v r e X r . 
The following stability result for it™ holds: 

Theorem 3.3. 27ie solution u" o/ ( 3. 12| ) satisfies the following bound 



(3.13) 



N-l 
n=Q 



/ 



n+l| 
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Proof. Choosing v r :— u" +1 in (3.12), we get 
1 



(3.14) 



At 



(u 



n+l 



<,< +1 ) + A« +1 ,< +1 ) = (/ 



By applying the Cauchy-Schwarz inequality on both sides of (3.14| and simplifying 
by ||< +1 ||, we get: 



(3.15) 



,rc+l| 



Kll < At||/ 



n+l I 



Summing from to N — 1 the inequality in (3.151, we get (3.13) 



□ 



In order to prove an estimate for \\u n — «™||, we will first consider the Ritz 
projection w r G X r of u G X: 



(3.16) 



A(u — w r , v r ) = Vty G X r 



The existence and uniqueness of w r follow from Lax-Milgram lemma. We now prove 
an estimate for u n — to", the error in the Ritz projection. 



Lemma 3.4. The Ritz projection w™ of u n satisfies the following error estimate: 

V 1 1 7/" -™ n ll < a\ (\ + a /ll Mr 1 IU + rv~ lN ) 



1 N 

(3-^) ^ E 



ra=l 



JV 



iV 



E ^ 

\ J=r+1 



\ 71=1 ^ j=r+l J ) 



Proof. Setting u :— u n in (3.16), we get: 

(3.18) A(u n - <, v r ) = Vu r G X r . 

We decompose the error u n — u>™ as u n — := (w n — /^(u™)) — (w™ — Ih tr (u n )) = 
f] n — </>", where Ih.r(u n ) is the interpolant of u n in the space X r . By the triangle 
inequality, we have: 



(3.19) 



1 N 

-Y\ 

N ^ 1 

n=l 



1 



N 



<n < ^ E ifo" 

71=1 



1 " 



We start by estimating ||?7 n ||. We note that Ih, r (u n ) consists of two parts: We first 
consider uJJ, the FE solution of (1.1), which yielded the ensemble of snapshots V 



defined in (2.1). Then, we interpolate u% in X r , which yields Ih,r{u n ). Note that 
this is different from [15] . where only the first part was present (see (8) in [15]). 



(3.20) 



1 N 

^Eii^i 



1 N 

- X> B --M« n )ll 



n=l 
1 N 

<-Y IK 

- N ^ 11 

71=1 



1 - 

^Ek-v(^")ii- 
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Using Assumption [2j it is easily shown [34] that: 

N 

(3.21) -^ik-<h<c^+ i -x: 



N 



n=l 



|m+l- 



Picking If l r (u n ) := V^(it n , ^jOx ¥>j m the last term on the RHS of (3.201 and then 
using (3.10), we get: 

N d 

(3.22) _J2\K-h,r(u n )\\< x Y, V 

n— 1 \ j=r-t-l 



Note that we consider that the time instances t n = n At in the time discretization 
(3.12 ) are the same as the time instances at which the snapshots were taken. If this 
is not the case, one should use a Taylor series approach (see (4.8) in [31]). 
Plugging ([3721]) and (|3~22"| in ([3720]) , we get: 



(3.23) 



JV N 

^Eii^i^ m+1 ^EiKiu+i + 



\ E ^ 



Similarly, using that X = Hq in (3.10), we get: 



(3.24) 



JV N 

^En^i^^E^i^ 



n=l 



\ E Ai 

\ ]=r+l 



Equation (3.18) implies: 

(3.25) A(u n - <, v r ) = A{ v n - #,« r ) = 0. 
Choosing v r — </>" in ( |3.25 1 implies: 

(3.26) A(#,#)=A(»A#f). 

We decompose the bilinear form A into its symmetric and skew-symmetric parts: 
A:= A s +A ss , where A s (u,v) := a {P' R Vu, P' R Vv)+s (Vu, Vv)+((g - |V • b) 
and A ss (u, v) :— (h ■ Vw + | (V • b) u,v) . Equation ( 3.26| ) implies: 



(3.27) A a {ft,W) + A 







A s ( v n ,W) + A ss (v n ,€) 



Assumption nl implies that A s ((/>™, </>") > C ||^>™||i )£ . Q - Thus, using the Cauchy- 
Schwarz and Young's inequalities, (3.27) becomes: 

c \m\l e , a < M€, #) 1/a Mv n , v n ) 1/2 + A ss ( v n , tf) 
(3 ' 28) < \ MW, €) + \ Mv n , v n ) + (b^,vo + \ ((V • b) r, n , o . 

Rearranging and using Assumption [I] ( |3.28 ) becomes: 
(3.29) C||CIIU« ^ C + VOI + I ((V • b) I 
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We now estimate each term on the RHS of (3.29) 
(3.30) 



|A a (»?W)|=e||WT + .g--V-b r, n ,r, 



1. 



a\\P R ^V n \\ 2 <C\\v n \\l £ , a 



nil 2 



To estimate the second term on the RHS of (3.29), we first note that ||P r]| < 1 
(since P R is L 2 -projection) and use the inverse estimate (3.5) in Lemma |3.l| to 
obtain: 

(3.31) iiiMvo < iivcii < n#iiiP < \fwFh\m\ ■ 

Using that {P R u, P R v) = Vu, v, the Cauchy-Schwarz and Young's inequalities, 
and the inverse estimate (3.9), we then get: 

(3.32) |(V\V#?)| < |(P ii (b 77 ™) ) Pfi(V^))| + |(P fl (b<),P fl (V^))| 

< \\P R (b V n )\\ ||P fl (V#)|| + \\P R (brf)\\ \\P R (Vr r )\\ 



< (^c\\M-%\\p R (bn n )f + ^\m\ 3 

J-||P>,f)f + f \\P R (VW)\\ 2 
2 a 2 

We note that this is exactly why we need the inverse estimate in in Lemma |3.1| 
to absorb ||<^|| 2 in the LHS of ( |3~29| ). If we had used ||V0; l || 2 instead, then we 
would have had to absorb it in e ||V<?!>"|| 2 on the LHS, and so the RHS would have 
depended on e. Finally, by using the Cauchy-Schwarz and Young's inequalities, the 
third term on the RHS of (3.29) can be estimated as follows: 

(3.33) | ((V • b) „« #) | < C WrW < C (~ ||r?"|| 2 + 1 1|<|| 2 



Collecting estimates ( |3.29[ ), ( |3.30[ ), ( |3.32[ ) and ( |3.33[ ), we get: 
(3.34) 



1 



ll^llLa < C[ h"||L, Q + - \\M-% \\P R (hr}* 



\v 



The last term on the RHS of ( [3~34] ), can be absorbed in C ||?7"|| 2 , £iQ . Since ||Pr|| < 1 
(P R is L 2 -projection) and ||b|| < 7 (by Assumption [lj, we get: 



(3.35) 



M-%\\P R (br, n )f<C\\M-%\\ V n \ 



n\\2 



Since ||Ph|| < 1 {P R is L 2 -projection) and ||b|| < 7 (by Assumption [T]) , we get: 



(3.36) 



2a a 
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Thus, using (3.35) and (3.36) in (3.34), we get 



(3.37) 



Wr\\\,e, a < C[h n \\ 2 + e || V^|| 2 + a HP^W 1 )!! 2 + C \\M-% II 



±\\P R (hr 1 n )\\ 2 + -\\ V n r 
Za a 



Since Pr is I/ 2 -projection, \\P R \\ < 1, and thus the second term on the RHS of 
(pH7| can be bounded as follows: a ||P fi (Vry")|| 2 < a ||V?/"|| 2 . Summing in (337}, 



wc get: 



(3.38) 



N N 
V E H^Hl,e,a < C (1 + \\M~% + a' 1 ) - ]T ||, ? "|| 2 



iV 



1 W 



n=l 



Using (3.23) and (3.24) in (3.38), we get 



1 N { 

(3.39) = ^ CHl + HM-ia + a- 1 

n=l l 



-l\l/2 



-li 

N 



N 

E 



U ||m+l 



\ n=l 



\ 



j=r+l 



IKIL+1 + 



, E ^ 

\ j=r+l 



Using (3.19), (3.231, and (3.39), we get (3.171 



□ 



Corollary 3.1. The Ritz projection of u n satisfies the following error estimate: 
(3.40) ||(tt"-<) t || < CUl + WM-iy + a- 1 ) 1 ' 2 



\ 



E >» 

j=r+l 



h m+1 \\u t \\ L2{Hm+1) + 
+Ve + a h m IKIUa^+i) + 



, E a. 

\ j=r+l 



Proof. The proof follows along the same lines as the proof of Lemma |3.4| with 
the error u n — replaced by (u n — u>") t = (r] n — 0™)t- Note that it is exactly 
at this point that we use the fact that the finite difference quotients du(t n ) are 
included in the set of snapshots (see also Remark 1 in [28]). Indeed, as in the 
proof of Lemma 3.4 the error (u n — w") t is split into two parts: (u n ~ := 



« - l h ,r(0) ~ (M" " hAO) = V? - As in (3.20)-(3.23), < can be 
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estimated as follows. 

\m<M-uu\+\\ui t -i h ,M)\\ 



(3.41) 



<C h m+1 \\u t \\ m+1 + 



, E 

\ j=r+l 



where in the last inequality in (3.411 we used (3.10) 



□ 



We are now ready to prove the main result of this section. 
Theorem 3.5. Assume that 

(3.42) L R = VX R = span{V(pi,...,V<p R }. 

Then the following error estimate holds: 

(3.43) 



N 



1 N { 

— £ \\u n - <|| < C (1 + \\M~% + a- 1 ) 172 

( 1 N 

U m+1 ^E(IK|| m+1 + |Kiu 2(Hm+1 )) + 

\ 71=1 

/ ! " 

+Ve + a \h m — (IKIL+i + \\uth*(H 
\ n=l 



\ E 



, E ^ 

\ J=r+1 



+||«°-«°|| +At||n t i|| L 2 (z ,2 ) 



\ 71=1 



|m+l 



v E 

\ j=ii+i 



Proof. We evaluate (3.11) at £ n +i, wc let i> = w r , and then we add and subtract 

ra+l _ u n 

a > v r 

At 



(3.44) 



,71+1 



At 



-,v, 



At 



+ a{u n+ \v r ) = {f n+ \v r ) 



Subtracting (3.12) from (3.44), wc obtain the error equation: 



(3.45) 



,71+1 



At 



-,v r - 



-,v, 



At '7 V At 
+ A(u n+1 - u? +1 ,v r ) + (a - A)(u n+1 ,v r ) = 0. 



We now decompose the error as u n — u™ = (u n — to") — (tt™ — u>") = ?y 
which, by the triangle inequality, implies: 



(3.46) 



IK -<||<lbf1 + ||^||. 



We note that has already been bounded in Lemma 3.4 Thus, in order to 

estimate the error, we only need to estimate The error equation (3.45) can 
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be written as: 



(3.47) 



,71+1 



At 



~~At ' 



/.n+l An 



+A( V n+1 - $ n +\v r ) + (a - A)(u n+ \v r ) = 



At 

n+l 



We pick v r := </>; l+1 in ( [317] , we note that, since </>™ +1 e X r , A(r? n+1 , <£™ +1 ) = 0, 
and we get: 

(3 4g) + ^ ~ = ~ (V n+1 ~ 

+ (r n ,ct>?+ 1 ) + (a-A)(u n+ \c%+ 1 ), 



We now start estimating all the terms in (3.48). 



where r n = u^ +1 - 

At 

The terms on the LHS of |3.48[ ) are estimated as follows: 

in+l m+h ^ /3IIJ.n+l||2 , _ || V7 j.n+1 II 2 , „, || D ' w j.ri+1 1|2 



(3.49) A(C +i , > P + e ||VC +1 || 2 + a ||P fl V# +1 H' 

(3.50) 



1 ^ n+1 -o; i+1 )> ~ (IIC +1 II 2 -IICIIIIC +1 II) 



Af 



Now we estimate the RHS of (3.48) by using the Cauchy-Schwarz and Young's 
inequalities: 



(3.51) 



1 

a7 



( V n+1 -77")+r™,^ +1 



< 



1 

a7 



(?7 n+1 - rf 1 ) + r" 



< 



2/3 



Af 



II0" +1 | 

m 



P n j.n+ll|2 



(3.52) 



(a - A)(u" +i , 0' r l+1 ) = -a (P R Vu n+i 1 P R V^ +i ) 
< a \\P R Vu^\\ \\P' R V^ +1 \\ < f ||P fl W l+1 || 2 + f ||P fl VC +1 f. 



Using ( 3.49 )-( 3.52) and absorbing RHS terms into LHS terms, (3.48) now reads: 
(3.53) 



^ (IIC +1 II 2 - iicn ||0; i+1 ||) + £ nc +1 || 2 + e ||v^ +1 n 2 



"ll^vc +1 ll 2 < ' 



2/3 



1 
At 



n+l _ ^n) + r r 



"\\P'M l+1 \\ 2 . 



By using Young's inequality, the first term on the LHS of (3.53) can be estimated 
as follows: 



(3.54) 



iic +i ii 2 - ncii W€ +1 \\ > ii0; i+1 n 2 - \ www 2 I W€ +1 w 2 
= \w T +1 t-\\m\ 2 - 



Using (3.54) in (3.53) and multiplying by 2 At, we get: 



(3.55) urT-Ur\\ 2 + ^t\\€ +1 \\i 



<C [At 



< C [At 



1 

~A~t 
1 

At 



in 



n+l _ r; «) + r r 



a At \\P R Vu 



Tt + l||2 



(rf +1 - v n : 



At\\r n \\ 2 + aAt\\P R Wu n+1 \\ 2 
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Summing from n = to n = N — 1 in (3.55), we get: 
N-l / 



N-l 



max IICH 2 + J2 A* IIC +1 HL Q < C At £ 



0<n<N 



(3.56) 



n=0 



At 



(rf + 1 - rf ) 



AT-l N-l \ 

U° r \\ 2 + At ]T \\r n f + a At ]T \\P' R Vu n+1 \\ 2 . 

n=0 n=0 / 



Proceeding as in [39] (see also [15]), we estimate the first term on the RHS of (3.56) 
as follows. We start by writing: 



(3.57) 



v n+1 -rT = 



t„+i 



Vt dt . 



Taking the L -norm in (3.57) and applying the Cauchy-Schwarz inequality, we get: 



\V 



n+1 „n I 



(3.58) 



< / l\\rit\\dt< 



l 2 dt 



1/2 



1/2 



< (At) 1 / 2 



t»+i \ V2 

2 



which implies At 



At 



n+l _ ^n- 



N-l 



N-l, we get At 



to n 

Corollary |3.1| We thus obtain 

N-l 



At 



< 



(v n+1 - v n ) 



1/2 

||7yt|| 2 dt I . Summing from n = 



(3.59) At ^ 



n=0 



1 

At 



(vf + 1 - rf) 



< ||//t||z,2(L2), which was bound in 

KC^l + WM^h+a- 1 ) 1 ' 2 
h m+1 \\u t \\ LHHm+1) + 

We + a k m ||ttt||joa(flm+i) + 



\ J=r+1 



\4 



r+l 



To estimate the third term on the RHS of (3.56), we use a Taylor series expansion 
of u n around u n+1 : 



(3.60) 



u n = u n+L - < +1 At + / u ft (s) (t„ - s) ds 



Taking the L -norm in (3.60) and applying the Cauchy-Schwarz inequality, we get 

rt»+i 



< 1 \\utt\\ ds < (At) 1/2 11^11^2(^2) . Summing from n = 0ton = N-l, 



we get: 
(3.61) 



N-l 



At J2 \\r n \\ 2 < At 2 \\u tt \ 



2 

L 2 (L 2 ) 



n=0 
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To estimate the last term on the RHS of (3.56), we use the fact that L R — S7X R 



(assumption ( 3.42 1 ) . We emphasize that this is the only instance in the proof where 
the assumption L H — V X R is used. Thus, we get: 



N-l 



N-l 



(3.62) a At ll^flVw" +1 || 2 = «At ^ ||Vu ,l+1 - P R Vu n+1 

n=0 n=0 

(pa) i "z* 

Ca- V inf ||Vu n+1 - Vv R \ 

( i N 



\ 3=R+1 



Using (3.59), (3.61), and (3.62) in (3.56), the obvious inequality max ||<jf>"|| > 

1 ' 1 ' ' ' 1 ' 0<n<N 



N 



1 N 

— T\ 

+ 1 i 

71=0 



inequality (3.46), and the estimates in Lemma 



error estimate (3.43) 



3.4 



we obtain the 
□ 



4. Numerical Results 



The goal of this section is twofold: (i) to show that the new VMS-POD model 



(2.9 ) is significantly more stable numerically than the standard POD-G model ( 2.7 1; 



and (ii) to illustrate numerically the theoretical error estimate (3.43). We also 



use Theorem |3.5| to provide theoretical guidance in choosing an optimal value for 
the artificial viscosity coefficient a and use this algorithm within our numerical 
framework. Finally, we show that the VMS-POD model (2.9) displays a relatively 
low sensitivity with respect to changes in the diffusion coefficient e. Thus, we 
provide numerical support for the theoretical estimate (3.431, which is uniform 
with respect to e. 

The mathematical model used for all the numerical tests in this section is the 



convection-dominated convection-diffusion equation ( 1.1 ) with the following param- 
eter choices: spatial domain fi = [0, 1] x [0, 1], time interval [0,T] = [0, 1], diffusion 
coefficient e — 1 x 10~ 4 , convection field b = [cos -|, sin ^] T , and reaction coefficient 
g = l. The forcing term / and initial condition uq(x) are chosen to satisfy the exact 
solution u(x : y,t) — 0.5 sin(7rx) sin(7r?/) [tanh ( ^^"o^ ) + l] > which is similar to 
that used in [13] . As in the theoretical developments in Section [3j in this section we 
employ the finite element method for spatial discretization and the backward Eulcr 
method for temporal discretization of all models investigated. All computations 
are carried out on a PC with 3.2 GHz Intel Xeon Quad-core processor. 



We start by comparing the VMS-POD model ( 2.9 ) to the standard POD-G model 



(2.7). To generate the POD basis, we first run a DNS with the following parameters: 



piecewise quadratic finite elements, uniform triangular mesh with mesh-size h = 
0.01, and time-step At = 10~ 4 . A mesh refinement study indicates that DNS mesh 

N 

E 



resolution is achieved. The average DNS error is 



l 

N+l 



= 2.04 x 10" 



71=0 



where TV = 1000, and u n and it^ are the exact solution and the finite element 
solution at t = nAt, respectively. The CPU time of the DNS is 9.42 x 10 4 s. Since 
the forcing term is time-dependent, the global load vectors are stored for later use 



VARIATIONAL MULTISCALE PROPER ORTHOGONAL DECOMPOSITION 15 



in all the ROMs. The POD modes are generated in H 1 by the method of snapshots; 
the rank of the data set is 104. For both POD-ROMs (POD-G and VMS-POD), 
we use the same number of POD basis functions: r — 40. 



Table 1 . Average errors for the POD-G model (2.7| with different 



values of r. Note that the POD-G model yields poor results. 



r 


20 


40 


60 


80 


N 

wnE h n -<\\ 

n=0 


1.25 x 10- 1 


1.11 x lO" 1 


9.28 x 10~ 2 


8.20 x 10~ 2 



We first test the POD-G model (2.7 1. The CPU time for the POD-G model is 



96.4 s, which is three orders of magnitude lower than that of a brute force DNS. 
The numerical solution at t = 1 is shown in Figure [I] for both the DNS (top) and 
the POD-G model (middle). It is clear from this figure that, although the first 40 
POD modes capture 99.99% of the system's kinetic energy, the POD-G model yields 
poor quality results and displays strong numerical oscillations. This is confirmed 

N 

by the POD-G model's high average error j^py — u r I = 1-H x 10 — 1 , where 

71=0 

it™ is the POD-G model's solution at t = nAt. Indeed, the POD-G model's average 
error is almost three orders of magnitude higher than the average error of the DNS. 
The average errors for different values of r listed in Table [I] show that increasing 
the number of POD modes (r) does not decrease significantly the average error. 
It is thus clear that the straightforward POD-G model, although computationally 
efficient, is highly inaccurate. 

Next, we investigate the VMS-POD model (2.9). We make the following pa- 
rameter choices: R — 20 and a = 4.29 x 10 . The motivation for this choice 



is given later in this section. The CPU time for the VMS-POD model (2.9) is 



106.2 s, which is close to the CPU time of the POD-G model ( |2.7[) . The numerical 
solution at t = 1 for the VMS-POD model is shown in Figure |1| (bottom). It is 
clear from this figure that the VMS-POD model is much more accurate than the 
POD-G model. Indeed, the VMS-POD model results are much closer to the DNS 
results than the POD-G model results, since the numerical oscillations displayed by 
the latter are dramatically decreased. This is confirmed by the VMS-POD model's 

N 

average error ^ — u™| = 4.48 x 10~ 3 , where u™ is the VMS-POD solution 

n=0 

at t = nAt: this error is more than 20 times lower than the error of the POD-G 



model. In conclusion, the VMS-POD model (2.9) dramatically decreases the error 



of the POD-G model (2.7) by adding numerical stabilization, while keeping the 



same level of computational efficiency. 

We now turn our attention to the second major goal of this section - the numer- 



ical illustration of the theoretical error estimate (3.43). Specifically, we investigate 



whether the asymptotic behavior of the RHS of estimate (3.43) with respect to R 



is reflected in the numerical results. We focus on the asymptotic behavior with 
respect to R since this is the main parameter introduced by the VMS formulation; 
the asymptotic behavior with respect to r was investigated in [7J, whereas the as- 
ymptotic behavior with respect to h and At is standard [HI [2S] ■ To investigate the 
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Figure 1. Numerical solution at t = 1: DNS (top), POD-G model 
( [27| (middle), and VMS-POD model ((2T9j> (bottom). Note that 
the VMS-POD model is much more accurate than the POD-G 
model, decreasing the unphysical oscillations of the latter. The 
CPU times for both the VMS-POD and POD-G models are three 
orders of magnitude lower than the CPU time for the DNS. 





d 

asymptotic behavior with respect to R, we have to ensure that ^/a, / ^2 <\j (the 

V j=R+i 

only term that depends on i?) dominates all the other terms on the RHS of ( |3.43 1 . 
To this end, we start collecting all the terms that depend on the exact solution u 



VARIATIONAL MULTISCALE PROPER ORTHOGONAL DECOMPOSITION 17 

Table 2. VMS-POD model's average error e = 

N 

Aril E — u r\\ an d its e 3 component for different values 

71=0 

of i?. 



i? 


e3 


1 ^ 
77+1 E M " 

n=0 


-<ll 


1 


1.29x10" 


i 


2.55x10" 


-2 


4 


9.34x10" 


2 


1.78x10" 


-2 


7 


6.69x10" 


2 


1.37x10" 


-2 


10 


4.68x10" 


2 


9.80x10" 


-3 


13 


3.20x10" 


2 


6.99x10" 


-3 



and we include them in the generic constant C. Next, we assume that the POD in- 
terpolation error in the initial condition — is negligible. We also assume that 
the time-step is small enough to neglect At \\uu \\l 2 (l 2 )- With these assumptions, 
the error estimate (3.43) can now be written as e < C (e\ + e2 + 63), where e is the 



VMS-POD model's average error, C a generic constant independent of r, i?, ft, At 
)h m+1 , 



and a, e\ 



\M: 



C'2 



IIm^ii; 




\j, and 



A,- 



d 

£ 

j=r+i y J=R+l 

To ensure that dominates the other terms, we choose r = 100 and consider 
relatively low values for R. This choice for r, which is not optimal for practical 
computations, ensures, however, that dominates e2- We also note that, when h 
is small, e% dominates e\ too. Thus, to investigate the asymptotic behavior with 
respect to R of the RHS of (13.431), we fix a — 5 x 10~ 3 , vary R from 1 to 14, 



and monitor the changes in e^. We restrict R to this parameter range to ensure 
that il ^2 (and thus e^) dominates ei and e\. Table 

V 3=-R+l 

JV 



lists the VMS-POD 



model's average error e 



N 



+1 ^„ 



it 



11,: 



and its e% component for different 



n=0 



values of R. We emphasize that, in this case, e% dominates the other two error 
components e x = 3.81 x 10~ 3 and &2 = 2-87 x 10 -3 . To see whether the theoretical 



linear dependency predicted by the theoretical error estimate (3.43) is recovered in 



the numerical results in Table [2] we utilize a linear regression analysis in Figure [2] 
This plot shows that the rate of convergence of e with respect to e$ is 0.9, which is 
close to the theoretical value of 1 predicted by (3.43). We believe that this slight 
discrepancy is due to the fact that the mesh-size h = 0.01 that we have employed 
in this numerical investigation is not small enough for our asymptotic study. 

Summarizing the results above, we conclude that the theoretical error estimate in 



(3.43 1 is recovered asymptotically (with respect to R) in our numerical experiments. 



Next, we use Theorem |3.5| to provide theoretical guidance in choosing an opti- 
mal value for the artificial viscosity coefficient a. The main challenge is that the 
theoretical error estimate (3.43) is asymptotic with respect to h, At and r, while in 



practical computations we are using small, yet non-negligible values for these pa- 
rameters. Furthermore, the generic constant C is problem-dependent and can play a 
significant role in practical computations. Notwithstanding these hurdles, we choose 
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Figure 2. Linear regression of VMS-POD model's average error 
with respect to e^. The convergence rate is 0.9, which is close to 



the theoretical value of 1 predicted by (3.43) 



10 



10 



0.02 



. P E = O D 1 V 5^° 



0.05 



0.1 



0.2 



a value for a that minimizes the RHS of (I3.43I): a = , H l =. In 



V j=r+l V i = R + 1 

the derivation of this formula, we made the same assumptions as those made in the 
numerical investigation of the asymptotic behavior of the VMS-POD model's error 
and we again considered that (3.43) can be written as e < C (ei + + e^). We 

/ d / d I d 

note that, if < / « 4/ an d h m « J tnen a becomes 

V ]=r+i y 3=ji+i y j=fi;+i 

too small in practical computations and the VMS-POD model becomes similar to 
the inaccurate POD-G model. To circumvent this, we use in our numerical tests a 
"clipping" procedure by setting a* = max {5, |}. 

Tablejsjlists the VMS-POD model's average error e = ^ £ ~ <ll for thc 

following values of r, R and a: r = 20, 40 and 60; R from 5 to r — 5 in increments of 
5; and a = 0.01a*, a*, and 100a*. Note that the VMS-POD model consistently 
performs best for a = a*. The only two slight deviations from this rule are for 
r = 60 (R = 20 and R = 30); we again believe that this is due to the mesh-size 



h = 0.01, which is not small enough for the asymptotic regime in Theorem 3.5 



Finally, we investigate numerically the VMS-POD model's sensitivity with re- 
spect to changes in the diffusion coefficient e. To this end, we run the VMS-POD 
model (2.9) with the same parameters as above (r = 40, R = 20 and a = a*) for 
different values of the diffusion coefficient: e — 10~ 2 , 10 -4 and 10~ 6 . Table [4] lists 
the average errors for DNS, POD-G and VMS-POD models for different values of 
e. It is clear from this table that the POD-G model's average error is significantly 
higher than the error of the DNS. The VMS-POD model, however, performs well 
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Table 3. VMS-POD model's average error e = 

JV 

]vTT S — u r\\ f° r different values of r and R, and 

71=0 

a = 0.01a*, a*, and 100a*. Note that the VMS-POD model 
consistently performs best for a = a*. 



r 


R 


0.01a* 




e 




a* 




e 




100a* 


e 






5 


1.2 x 10" 


:', 


1.0 x 10" 


l 


1.2 x 10" 


i 


5.8 x 10" 


■2 


1.2 x 10 1 


7.8 x 10" 


2 


20 


10 


2.0 x 10" 


-3 


9.5 x 10" 


2 


2.0 x 10" 


i 


2.4 x 10" 


-2 


2.04 x 10 1 


2.6 x 10" 


2 




15 


3.3 x 10" 


-3 


8.2 x 10" 


2 


3.3 x 10" 


i 


2.0 x 10" 


-2 


3.3 x 10 1 


2.5 x 10" 


2 




5 


6.4 x 10" 


; > 


1.09 x 10 


-1 


6.4 x 10" 


s 


3.0 x 10" 


2 


6.4 x 10 _i 


7.2 x 10" 


■2 




10 


1.1 x 10" 


-4 


1.0 x 10" 


1 


1.1 x 10" 


-2 


1.8 x 10" 


-2 


1.1 x 10° 


2.5 x 10" 


2 


40 


20 


4.2 x 10" 


-4 


9.7 x 10" 


2 


4.2 x 10" 


-2 


4.4 x 10" 


-3 


4.2 x 10° 


4.1 x 10" 


3 




30 


1.7 x 10" 


-3 


6.8 x 10" 


2 


1.7 x 10" 


1 


8.1 x 10" 


-3 


1.7 x 10 1 


1.0 x 10" 


2 




35 


3.0 x 10" 


-3 


4.9 x 10" 


2 


3.0 x 10" 


1 


2.1 x 10" 


-2 


3.0 x 10 1 


2.4 x 10" 


2 




5 


5.0 x 10" 


5 


8.7 x 10" 


2 


5.0 x 10" 


:i 


1.8 x 10" 


2 


5.0 x lO" 1 


7.0 x 10" 


2 


60 


10 


5.0 x 10" 


-5 


8.7 x 10" 


2 


5.0 x 10" 


-3 


1.3 x 10" 


-2 


5.0 x lO" 1 


2.4 x 10" 


2 




20 


5.0 x 10" 


-5 


8.7 x 10" 


2 


5.0 x 10" 


-3 


1.0 x 10" 


-2 


5.0 x lO" 1 


3.9 x 10" 


3 




30 


1.2 x 10" 


4 


8.0 x 10" 


2 


1.2 x 10" 


2 


4.4 x 10" 


-3 


1.2 x 10° 


7.4 x 10" 


4 




40 


5.4 x 10" 


4 


5.5 x 10" 


2 


5.4 x 10" 


2 


1.2 x 10" 


-3 


5.4 x 10° 


2.4 x 10" 


3 




50 


1.8 x 10" 


-3 


2.6 x 10" 


2 


1.8 x 10" 


1 


1.3 x 10" 


-2 


1.8 x 10 1 


1.4 x 10" 


2 




55 


2.9 x 10" 


-3 


2.2 x 10" 


2 


2.9 x 10" 


1 


1.1 x 10" 


-2 


2.9 x 10 1 


1.2 x 10" 


2 



for all values of e and displays a low sensitivity with respect to changes in the dif- 
fusion coefficient. Thus, we provide numerical support for the theoretical estimate 



(3.43), which is uniform with respect to e. 



Table 4. Average errors of DNS, POD-G and VMS-POD models 
for different values of the diffusion coefficient e. The POD-G model 
performs poorly. The VMS-POD model performs well and displays 
low sensitivity with respect to changes in e. 



£ 


DNS 


POD-G 


VMS-POD 


N 

7&i £ \\ un - u h\\ 


N 

ivtiE IK-<II 


]vTTEIh n -<ll 


lO -2 
10" 4 
10" 6 
10" 8 


1.10 x 10~ 4 
2.04 x 10~ 4 
1.88 x 10~ 4 
2.46 x 10~ 4 


1.10 x 10^ 

1.11 x 10" 1 
1.17 x 10" 1 
1.17 x 10 _1 


4.05 x W-" 4.27 x lO" 3 
4.29 x 10" 2 4.48 x 10" 3 
9.65 x 10 -2 4.05 x 10 -3 
1.01 x 10 _1 4.05 x 10 -3 



5. Conclusions 

We presented a new VMS closure modeling strategy for the numerical stabi- 
lization of POD-ROMs of convection-dominated equations. The new POD-ROM, 
denoted VMS-POD, utilizes an artificial viscosity term to add numerical stabiliza- 
tion to the model. Following the guiding principle of the VMS methodology, we 
only add artificial viscosity to the small resolved scales. Thus, no artificial viscosity 
is used for the large resolved scales. The POD setting represents an ideal framework 
for the VMS approach, since the POD modes are listed in descending order of their 
kinetic energy content. 
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A thorough numerical analysis for the finite element discretization of the new 
VMS-POD model was presented. The numerical tests showed the increased nu- 
merical stability of the new VMS-POD model and illustrated the theoretical error 
estimates. We also employed the theoretical error estimates to provide guidance in 
choosing the artificial viscosity coefficient in practical computations. We emphasize 
that the theoretical error estimates were uniform with respect to e, the diffusion 
coefficient. The numerical tests confirmed the theoretical results: The average error 
of the VMS-POD model showed a low sensitivity with respect to changes in e. 

Although the new VMS-POD model targets general convection-domainted prob- 
lems, it was analyzed theoretically and tested numerically by using the convection- 
dominated convection-diffusion equations. We chose this simplified mathematical 
and numerical setting as a first step in a thorough investigation of the new VMS- 
POD model. Next, we will utilize the new VMS-POD model in the numerical 
simulation of turbulent flows, such as 3D flow past a circular cylinder [40j . We 
also note that, to our knowledge, this is the first time that the VMS formulation 
used in |29j for the numerical stabilization of finite element discretizations has been 
used in a POD setting. We will investigate in a future study the alternative VMS 
formulation proposed in (T3] and compare it with the VMS-POD model that we 
introduced in this report. 
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